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A  - 


SIGNIFICANCE  AND  XP  LA  NATION 


Many  uthemafica  models  for  heat  flow  or  fluid  flow  involve  the  specifi¬ 
cation  of  a  flov.  i  ite  act  .s  the  boundary  of  a  region  which  may  depend  in  a 
i  tr  fashion  upon  the  unknown  variable  (e.g.  temperature).  Formulation 
and  anal  sis  of  efficient  numerical  procedures  for  approximating  the  solutions 
of  such  problems  are  studied.  Previously,  finite  element  methods  used  for 

modeling  these  physical  problems  have  boun  at  most  second  order  correct  in  the 

_i  .  it  4 


time-discretization  error.  Wo--f»rndiKw  methods  Jhich  are  second,  third,  and 

T  A 

fourth  order  correct  in  time  and  which  convert  the  nonlinear  problems  into 

solution  of  large  systems  of  linear  equations  via  an  extremely  stable  algorithm 

with  essentially  no  restrictions  between  sizes  of  time  and  space  discretizations. 

The  basic  multistep  methods  presented  produce  different  systems  of  linear 

equations  at  each  time  step.  A  preconditioned  iterative  stabilization  procedure 

is  presented  and  analyzed  which  allows  for  the  factorization  of  only  one  large 

matrix  to  be  used  at  each  time  level  in  the  solution  process.  Optimal  order 

(v 

error  estimates  are  obtained.  The  paper  also  contains  work  estimates  which 
show  the  large  computational  savings  of  the  preconditioned  iterative  stabiliza¬ 
tion  technique.  Almost  optimal  order  work  estimates  are  obtained.  _ 


) 
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KI-T'I i '  1 ENT  MULTISTEF  PROCEDURES  FOR  NONLINEAR  PARABOLIC 
PROBliEMS  WITH  NONLINEAR  NEUMANN 
BOUNDARY  CONDITIONS 


Richard  E.  Ew i ng ' 


I.  Introduction.  We  shall  consider  the  numerical  solution  of  nonlinear  para- 
bol ic  partial  differential  equations  with  nonlinear  Neumann  boundary  conditions  of 
the  form 


a) 

c(x,u) 

3u 

3t 

V  •  [a (x. 

u) Vu  +  b(x,u)l  *  f(x,t,u),  x 

t  12,  t  £  J  , 

b) 

a(x,u) 

3u 

3v  + 

b(x,u)  • 

v  *  g (x , t  ,u)  ,  x  e  3S1,  t  c 

J  , 

(1.1) 

c) 

u  (x  ,0) 

uo 

(x)  ,  X 

.  ft  , 

where  ft 

is  a  bounded 

domain  in 

!Rl*  ,  d  •'  3,  with  boundary 

3ft,  v  is  the 

outward 

unit  normal  to  3ft 

I,  J 

( 0  ,  T 1  , 

and  c ,  a ,  b ,  f ,  g ,  and  u^ 

are  prescribed. 

We 

■■hall  use  a  Galerkin  approximation  in  the  space  variable  and  high-order,  efficient, 
multistep  time-stepping  procedures.  We  first  present  basic  multistep  time-stepping 
procedures  which  produce  a  different  linear  system  of  equations  to  be  solved  at 
each  time  step.  We  then  modify  the  basic  procedures  by  using  a  preconditioned 
iterative  method  to  approximate  the  solution  of  the  linear  equations.  The  use  of 
a  t ime- independent  preconditioning  matrix  eliminates  the  need  to  refactor  a  new 
matrix  at  each  time  step,  while  the  iterative  procedure  stabilizes  the  resulting 
algorithm.  Using  this  modification,  we  obtain  the  same  order  error  estimates  as 
for  the  base  scheme  with  greatly  reduced  computational  requirements.  We  obtain 
very  nearly  optimal  possible  work  estimates  for  our  procedure. 

Galerkin  procedures  for  parabolic  problems  with  nonlinear  Neumann  boundary 
conditions  were  first  considered  by  Douglas  and  Dujiont  in  l HI  .  Then,  in  [171, 
Luskin  extended  this  work  of  [ H 1  to  quasilinear  equations  similar  to  those  con¬ 
sidered  here.  Luskin  used  Crank-Nicolson  time-stepping  methods  which  are  second 
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order  correct  in  the  time  discretisation.  In  11.']  ,  t  li«-  author  used  the  iterative 
stah  1 1  tzat  ion  techniques  developed  in  |'>,  1  01  to  present  compu  t  .t  t  i  ona  1 1  y  et/i.  lent 
variants  of  the  methods  of  Luskin  and  extended  these  methods  to  treat  coupled  ■  y1- 
terns  of  nonlinear  part dal  differential  equations  with  nonlineai  Uuindary  cond i t  ions . 
In  this  paper,  we  present  t ime-steppinq  procedures  which  are  hioher-otdet  m  time 
than  those  analyzed  in  lb-ll,  l  ],  These  t ime-st eppino  schemes  ate  based  on  the 
backward  differentiation  multistep  schemes  Jef.  IS,  14,  l'>| .  They  have  been  pre¬ 
sented  and  analyzed  for  quasilinear  parabolic  equations  by  bramble  and  Sammon  in 
12,  7].  Very  efficient  alternating  direction  variants  for  use  on  roetanqulai 
domains  will  appear  in  14,  S I . 

The  efficient  time-steppinq  techniques  presented  here  can  also  bo  used  to 
analyze  approximation  procedures  for  initial  boundary  value  problems  for  many  othet 
types  of  nonlinear  partial  differential  equations.  The  author  has  applied  itera¬ 
tive  stabilization  techniques  to  equations  of  Sobolev  type  (in  J 1 0 } >  which  have 
applications  in  thermodynamics ,  fluid  flow  in  fissured  rock,  and  shear inu  of  second 
order  fluids.  In  111,  12],  the  methods  are  applied  to  coupled  systems  of  equations 
which  model  miscible  displacement  in  porous  media.  Also  the  author  has  used  iteiu- 
tive  methods  successfully  for  second  order  in  time  equations  (in  ll'l)  which  have 
applications  in  vibrational  problems  and  nonlinear  viscoelasticity. 


In  Section  2  we  introduce  certain  notational  preliminaries  and  present  the 
base  time-steppinq  Galerkin  schemes.  In  Section  t  we  present  om  iterative  modi¬ 
fications  of  the  base  methods  and  analyze  the  effect  of  the  iterative  approximation 
on  a  sinqle  time  step.  In  Section  4  we  obtain  qlobal  error  estimates  for  a  partic¬ 
ular  multistep  method.  Section  5  contains  a  brief  discussion  of  the  computational 
complexity  of  the  methods  presented. 


II.  Preliminaries  and  Description  of  dalerkin  Methods,  bet  ( <f,  ijd  *  /  v'ydx , 

II ^11*  -  (*,iM  -  /3ij^i|ids,  an 

space  on  il  with  norm 


(2.U 


and  |t|>|  -  |*| 

W2 


|l|»|  "  0/’ 

,ip).  Let  W 

y  IN 

\  1/8 

1 1  <k|l 3xn  H 

8  (17)/ 

<*.  When 

s  -  2,  let 

(Fj , F^)  , 

write  ||  VF  II 

Ilf,  I!8.  ♦  ||  F  II* 

1  wk  2  wk 

s  s> 


1/s 


.  For  definitions  of  correspondinq  fractional  order  spaces, 


see  { 16] . 

Let  (m^)  be  a  family  of  finite-dimensional  subspaces  of  ll' (ill  with  the 
followinq  property: 

For  p  ■  2  or  p  •  ",  there  exist  an  inteqor  r  '  and  a  constant  K 


such  that,  for  1  <  q  <  r  and  <!’  t  w  (>J)  , 
-  P 


inf  {H*-\ll  +  hiu-xlt  }  <  K  11*11  hq  . 

X- M  w''  W*  °  w'1 

h  p  p  p 


.\V  also  assume  that  (til  satisfies  the  following  so-called  "inverse  assumption.  " 


,f  * .  v 


a)  H*H  ■  h  K  ll*H 
1  —  0 


b)  !*|  <_  h  "kqII *H  , 


-)  11*11  ♦  hllv*ll  <  k  h  "’ll*!! 

a>  V*3  —  (1 

L  <«)  L  (.:) 


Restrict  il  as  follows  (with  (?)  denoting  the  collection  of  restrictions): 
2 

1)  P.  is  H  -roqular. 

(?)  :  2)  3i'.  is  l.ipschitz. 

3)  There  exists  a  constant  such  that 

|*|2  <  KM  .  4  > 


If  X  is  a  normcd  space  on  Q  with  norm  11*11  and  *  :  1 0 , T 1  *  X,  then  we 

■fine 


a)  Ml 


L  ( J ;  X)  1(1 


/THv*(t) 


1 1/s 

'xdt 


1  <  s  <  “’  , 


b)  M 


sup  II*  (t)  I 


L,  ( J ;  X )  t.  10, T) 


Throughout  the  paper  we  shall  assume  that  a  and  c  are  bounded  above  \nd 
below  by  positive  constants  and  that  a,  b,  c,  and  g  are  smooth  functions  ot 
their  arciuments.  We  shall  also  assume  that  the  solution  u  is  sufficiently 
smooth  for  our  arguments  to  hold.  For  typical  explicit  smoothness  assumptions  on 
u  ar.d  the  coefficients,  see  [8-12,  17). 

As  iij  [18],  we  shall  introduce  an  auxiliary  elliptic  problem  to  aid  in  oui 
•  nil. sis.  I,et  \  >  0  be  chosen  sufficiently  larqe  that  the  bilinear  form 

N  (*.V,  X)  (a  (*)  V*  ,  VX)  +  \(*,X)  -  <q  (t  ,** )  ,X) 


:«i  r  i  ;  I  ics 


N  (*;*,* )  -  K0M1‘,  *,*  f 


Let  W  «  be  the  projection  of  u  into  ft  ,  defined,  for  each  t  e  J,  by 


N(uC.t);  WC,t),X)  -  N(u  ( • ,  t)  j  u  ( • ,  t )  ,X) 


(2.6) 


-  (c(u)  !pX>  ♦  (b(u),Vx)  ♦  ( f  (u)  ,  X)  ♦  Mu.X),  X  *  ^ 


Then,  as  in  18,  9,  12,  181,  we  can  obtain  the  following  lemma. 


Lemma  2.1.  There  exists  a  constant  (u)  such  that  if  n  ■  u  -  W,  s  =  0 

or  s  ■  1 ,  and  2  <_  q  r , 

a)  II  nil 


<  K  hq~*ll  ull 

®  s  —  1  m  a 

L  < J; H  )  L  (J»Hh) 


(2.7) 


b> 


3n 


-  „  ,  <  K,  hq_S{ll  ull  , 

lt  L2(J;HS)  1  L2(J;Hq) 


•  1 

3u 


In  order  to  require  weak  smoothness  assumptions  on  we  shall  need  to  use 

some  duality  theory  and  obtain  some  approx imat ion  theory  results  in  negative- 
indexed  norms.  For  these  results,  assume  that  ft,  a,  b,  c,  and  q  are  suffi¬ 
ciently  smooth  1161  that  for  each  t  <  J,  if 


a) 

b) 


la(x.u)Vu)  ♦  X^u  ■  <1^,  x  e  $5  , 


(2.8) 


/  ,  3u 

a(x.u)  ^  -  *2. 


x  f  , 


then 


,U^2  -  K(U>{|t*lllk  +  l^l  1}  ' 


(2.9) 


k+2 


V  4  2 

If  (2. 8) -(2. 9)  holds,  v#e  shall  say  that  (1  is  K  -regular.  Next,  define  for 
k  >  0, 


a)  ll*ll_k  =  supU*,*)  s  |l^llk  -  11  , 

b)  M  =  sup{<*,*>  :  |*|  -  1)  . 


(2.10) 


Lemma  2.2.  If  (1  is  H**2-reqular  for  k  <  1,  there  exists  a  constant  K(u) 
such  that  for  1  <  q  <  r  and  t  <  J, 


II  nil 


-k  *  |n|  1  4  |l?|  iK(u)hq+NHullq  4  |U|  >  .  (2.11- 

- 1  st  y)  -X  q 


Proof ;  See  112) 


We  also  make  the  assumption  on  (M.  1  and  u  that  there  exists  a  const  a: 


ant 


K,  such  that 


-4- 


oo  cv 

1-  ( J ;  L  ) 


II  vw  II 


a>  a> 

L  (J»L  ) 


3w 

II  3t 


*  "vl? 


L  (J;L  ) 


L  ( J ;  L  ) 


r’w 

2  I 

3t  n  “  1 

L  (J;ir) 


3  'w 
3t3 


L  (J;H  (3I!|) 


34W 

i 

3t' 


II  »  1 

L  ( J  ;  H  ) 


iK2 


'ufficient  conditions  for  the  above  to  hold  can  be  found  in  19,  10,  18]. 


We  next  consider  discrete-time  Galerkin  approximations.  Let  At  >  0, 

„  ..."  ,  .,.n,_  .  ...  .n, 


T/At 

<  X  and 

0 

t  =  oAt 

a) 

,  .  n+1 
d 

,  n+ 1  .  n 

i},  _,j, 

At 

b) 

. .  n+1 

<5tp  = 

,  n+1  n 

*l>  -  U> 

c) 

,2  n+1 

5  i|)  « 

,  n+1 

-  2* 

d) 

.3  n+1 

5  ip  = 

, n+ 1  ,  .  i 

e) 

c4  n+1 

6  \p  * 

,n+l  .,1 

ip  -  4\p 

n  ,  n-1 

+  t), 

n  ,  ,n-l  ,n-2 
n  .  , , n-1  , , n-2  .  , n-3 


(2.12) 


We  next  define  a  family  of  extrapolated  coefficient  backwards  differentiation  multi- 
step  discrete  time  methods. 

Let  U  -.  {tQ,...,t  }  -*  ft  be  an  approximate  solution  of  (1.1).  Assume  that 
U  are  known  for  k  <_  n.  Then,  given  certain  choices  of  parameters  0,  a^  ,  a,, 
a^,  and  a^  and  an  extrapolation  Un+^,  we  determine  Un+*  to  satisfy 

(c(0n+1)  U  -- t~U  ,  X)  +  8(a(un+1)VUn+1 ,  VX) 

=  6<g(tn+1,  Un+1),  x >  +  (c(Un+1)  j-  Ia.un+a0un"1+a,un'2+a.un_2l ,  X)  (2.13) 

At  1  2  3  4 

-  6(b(un+1),  VX)  +  8(f(tn+1,  Un+1),  X),  x  t  \  • 

A  particular  example  from  this  family  of  methods  is  the  choice  Un+^  *  un,  P  =  1 
and  =  0,  i  »  1,2, 3, 4.  This  choice  is  the  the  well-known  backward  Euler  method 

with  laqged  coefficients  which  is  known  to  have  time-discretization  error  of  order 
At.  Other  choices  of  the  parameters  and  extrapolation  in  the  coefficients  yield 

2  3  4 

temporal  errors  of  order  (At)  ,  (At)  ,  and  (At)  . 
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We  present  these  special  choices  in  the  following  table. 


Table  1 :  Selected  Multistep  Methods 


Extrapolat ion 

e 

*1 

a2 

°3 

U4  S 

" 

. . Jj, 

Time-discret izat ion 

Error 
.  (AU 11 

UnU-6UIl>1 

i 

0 

0 

0 

| 

0  | 
u 

At 

Un+1-62Un+1 

2/3 

1/3 

-1/3 

0 

0 

(At)2 

u“*l-«V*1  ~ 

6/11 

7/11 

-9/11 

2/11 

0 

fl 

(At)3 

Unfl-{4Un+1 

12/25 

23/25 

-36/25 

16/25 

- r 

-3/25  | 

(At)4 

We  note  that  by  extrapolat ing  the  coefficients  in  (2.13),  we  have  reduced  each  of 
the  above  problems  to  the  solution  of  a  different  set  of  linear  equations  at  each 
time  step. 

III.  Iterative  Stabilization  Procedures.  In  this  section  we  consider  effi¬ 
cient  methods  for  solving  the  linear  eauations  arising  from  (2.13).  We  note  that 
the  coefficient  matrices  from  (2.13)  change  with  each  time  step.  In  order  to  avoid 
the  factorization  of  different  matrices  at  each  time  step  to  solve  the  different 
systems  of  linear  equations,  we  shall  present  an  iterative  method  for  approximating 
their  solution  to  sufficient  accuracy. 

Let  be  a  basis  for  M^  and  let  um  from  (2.13)  be  written  as 

Um  -  T  e"V  .  (3.1) 

i-1 

Using  (3.1),  (2.13)  can  be  written  as 

4 

b  (0(5  “5  )  =  1C  (0{  l  a. 5  )  ♦  AtF  (5)1  (3.2) 

i=l 

H  Fn(0 

where  the  matrices  and  vectors  are  of  the  form 


a)  L  (•'.) 


»n  ,  n 
*  At  PA 


b)  Cn( :.)  -  (  (c  (  "  £+.)*.,*.))  . 

kil  k  k  j  1 


c)  An(t.)  -  ((a(  V  )V;.  ,Vv».))  , 


k  1  k  1  1 


a)  F! -  P(<q(tn+1,  V  £%  ),*' .)  -  (b{  l  £  V  ), 

k-1  1  k-1 


♦  (f(tn+1.  I  l'' *),*.))  , 

k-1  "  K  ‘ 


for  i,  j 


1 . m. 


Instead  of  solving  (3.2)  exactly,  wo  shall  approximate  its  solution  by  usiw 

0 

la  iti'ia*  ivr  procedure  which  has  boon  i  reconditioned  by  I.  ,  t  he  associated  mat  i  i>. 
with  coefficients  evaluated  at  t  =  0,  for  each  time  step.  The  preconditioning 
orocess  eliminates  the  need  for  factoring  now  matrices  at  each  time  stop,  while  the 
i*  -rative  procedure  stabilizes  the  resulting  problem.  The  stabilization  process  re- 
uires  iteration  only  until  a  predetermined  norm  reduction  is  achieved. 

IV no to  by 

v*  •  l 

k-1  k 

the  ipprox  imat i on  to  1'  produced  by  only  approximately  solving  (3.2).  An  itera¬ 
tive  procedure  for  obtaining  the  necessary  Vk  starting  values  using  the  iterative 
procedure  described  here  will  appear  in  13].  We  assume  such  a  starting  procedure 
i as  been  used  to  obtain  sufficiently  accurate  (see  (4.7))  starting  values.  Thus 

assume  V  ,  ...,V’n  have  been  determined.  We  shall  determine  the  M-dimensional 

vector  y!‘  ‘  (and  thus  V°+^)  using  a  preconditioned  iterative  method  to  approx i- 

n+1  n+1  n 

rate  from  (3.2).  As  an  initial  guess  for  £  -  t.  we  shall  extrapolate 

tram  previously  determined  values.  Specifically,  for  a  particular  method  having 

t.ime-t : uncation  error  (At)',  we  shall  use  as  the  initialization  for  our  iterative 

procedure 

X  =  (Y"+1  .  5»+V+1  ,  (3.5) 


'cio  the  backward  difference  operator  i5m  is  defined  in  (2.12)  for 

i ,  .  .  . ,  4  .  Since  we  are  usinq  previously  determined  > 1  in  the  coefficient 
^  ■  n+1  _  .  _  ,  _ . 


•ices  to  determine  > 


our  errors  accumulate 
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In  order  to  estimate  the  cumulative  error,  we  first  consider  the  single  step 
-n+1 


error.  We  define  y 


to  satisfy 

Ln(»(Yn+1  -  yn)  *  Fn(Y),  n  >  V 


(3.6) 


We  can  use  any  preconditioned  iterative  method  which  yields  norm  reductions  of  the 
form 


i  n,  .1/2  -n+1  n+1,,.  _  n.  .1/2  -n+1  n+1  y+1  n+1.,, 

lL  (Y)  (Y  -  Y  )H  <  P  II L  (Y)  (Y  -  Y  +  6  Y  )  II 

e  -  n  e 


(3.7) 


where  0  <  <  1  and  the  subscript  e  denotes  the  Euclidean  norm  of  the  vector. 

A  specific  iterative  procedure  for  obtaining  (3.7)  is  the  preconditioned  conjugate 
gradient  method  analyzed  in  11,  9,  10). 


Let 

a)  IMI2  =  (c(Vn+1)se,*>)  , 

n 

c 

b)  IMI2  =  (a(Vn+1)V*,W!)  , 

n 

a 

O  IIMII  =  IMI  +  (At)1/2IMI 

'  n  n  i 

c  a 


be  special  norms  and  seminorms.  Note  that  II  •  II  and 

n 

c 

valent  to  11*11  and  117*11,  respectively.  Then  letting 


(3.8) 


II  are  uniformly  equi- 
a 


M 


v*  .  y  v"*  , 

L  '  i  i 
i=l 


(3.9) 


with  y*"  defined  in  (3.6),  we  see  that  v1141  satisfies 


(c(vn+1)  V-  — ■  -,X)  +  8(a(Vn+1)VVn+1,VX)  +  6(b(Vn+1) ,VX) 
At 


(3.10) 


=  6<g(tn+1,vn+1) ,X>  +  B(f (tn+1,Vn+1) ,X)  +  (c(Vn+1)  l  a.Vn+1-1,x),  X  e  ^  • 

i=l 


Also  using  (3.8),  our  single-step  error  (3.7)  becomes 

t||vn*'1-v"+l|||  <  Jj!_  |||5»*V*1|||  n>„*l  .  (3. 

Mn 

We  note  that  as  in  16,  12) ,  there  is  a  Q  depending  upon  bounds  for  the  coeffi¬ 
cients,  such  that 
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r 


cith  <  Q  •  1,  ami 

nAt ,  n  ■  1 


(3.12) 


IV. 


A  Priori  Error  Estimates.  In  this  section  we  develop  a  priori  bounds  for 


• h>  errors  V*‘  -  u'  for  the  procedures  defined  in  (3.10)  using  the  base  schemes  de- 
r mce  in  (2.13).  The  techniques  for  treating  the  nonlinearities  in  the  coefficients 
>t  i,  b,  and  f  are  tedious  and  appear  in  (7,  9,  12].  Therefore,  for  simplicity 
•t  exposition,  we  shall  consider  the  simplified  problem 


a)  clx,t) 


3u 

3t 


!  a  ( x  ,  t )  V  u )  =  0 


3u 

b)  a (x , t)  —  =  q(x,t,u) 

3v 


c)  u(x , 0)  =  u  (x) 

0 


X  (  .  ,  t  C  J 

X  £  •■  ,  t  €  J 

X  £  0 


(4.1) 


We  an  thus  examine  the  hioher-order  efficient  time-stepping  procedures  without  the 
adrit  !  complexity  of  nonlinearities,  except  in  the  Neumann  boundary  condition. 

Also,  for  simplicity,  we  shall  present  the  details  for  the  particular  method 

pi 

whose  choice  of  parameters  yields  time-discret ization  error  of  order  (At)  where 
..  =  3.  Proofs  of  stability  and  convergence  for  the  other  methods  follow  similarly 
and  can  be  derived  from  the  proofs  of  similar  problems  which  appear  in  17] . 

For  v  =  3 ,  the  base  approximation  scheme  for  (4.1)  from  (2.13)  can  be  writ¬ 
ten  as 


(c  «un+1,X)  +  ~  At (a  ?un+1,VX) 
n+1  11  n+1 


At(g(tn+1,un+1)  ,X>  +  <cn+ll~  Sun-  yy  {Un"Al  ,X)  ,  Xe 


n-1. 


where  c 


n  +  1 


/  .n+1. 

c  (x, t  ) ,  a 


n+1 


,  .n+1.  j  .'.n+1  ?T,n  . ,n — i  ,  ,,n-2 

a(x,t  )  ,  and  U  =  3U  -  3U  +  U 


(4.2) 


Let 


n!'  =  un  -  wn  and  £n  =  Vn  -  wn.  We  know  from  Lemma  2.1  that  W  is  a  function  in 
which  is  sufficiently  close  to  u.  We  next  estimate  how  close  V  and  W  are. 

From  (2.6),  (4.1),  and  (4.2),  we  obtain  the  followina  error  equation 


-g- 


'Vl'4"*1-'1  ♦  TT(*n+lV;n+1'Vx) 


,n-l 


n+1 


11 


«C  1  , X)  +  lYY  AAt(n"T\x) 


+  At (c  Id  n 
n+1  t 


n+l  7  n  2  n-1,  .. 

•  rrdtn  *  it  V  '•xl1 


n+1 


-  6un+1  +  £  6un 


2  .  n-1.  . 

YY  6u  1 'X) 


*£-4t<»(t"*l,Vn*1>  -  ,(t"+1.»n+1).x> 


*  l<cn*l(V"*'  '  v"*1’ ’V  *  n-  st(«n+iv<vn*1  -  v"*') .5x1] 
-  /  ,  7  ,„n  2  -,n-l,  .  n+1  n+1  ,  , 

Cn+1  Tl  6?  '  IT  ^  lrX)  +  T1  lX)  +  T2  (X) 


+  t”+1(x)  +  t"+1(x),  X  e  Mb 


(4.3) 


Tern  enters  because  we  are  comparing  V  to  W  instead  of  directly  to  u. 

Term  T  measures  how  well  the  multistep  scheme  approximates  —  and  term  T 

2  ,  dt  3 

arises  from  the  nonlinearity  of  g.  Finally,  the  single-step  error  made  by  using 
the  iterative  procedure  to  approximately  solve  the  linear  equations  appears  in  term 
T, . 


We  shall  first  present  a  few  lemmas  which  will  help  separate  the  various  parts 
of  our  analysis.  First  we  note  that  the  parameters  8 (p )  and  a^(p),  i  =  1,...,4, 

are  chosen  in  (2.13)  to  insure  the  following  consistency  result. 

Lemma  4.1.  For  each  p  =  1,2, 3, 4,  the  choice  of  parameters  0(p)  and  a. (p) , 
i  =  1,2, 3, 4  given  in  Table  1  yields 

n+1  ,4 

II  0  (p)  At  -  I6un+1  -  l  a.(y)u ]  II  <  K3(At)P+1  .  (4.4) 

i=l 

We  next  consider  the  following  lemma  which  will  provide  the  estimates  for  the 
basic  stability  of  our  methods. 

Lemma  4.2.  Assume  that  Z°  satisfies,  for  m  >_  2 , 

J.-1 

I  I(cn+16zn  ,X)  +  yj-  At(an+1VZn+1,7X)] 
n=m  (4.5) 

£-1 

»  l  f(cn+1t^  azn  -  Yi  6zn_1l'x>  +  <Fn+1rX)i,  X  €  Mb  . 

n=m 
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hen  there  exist  constants  K.,  K_  and  K  such  that  setting  X  =  zn+^  yields 

4  5  fa 


\?.if  +  Yi«^n+1ir  +  izn+iif  Ati 

n+l 

n=m  a 


(4.6) 


-  K.  [1(2  II  4 

—  4 


Y  iu2n+i 

n-m-2 


f-1 


+  l  llzn+1||At  +  |  l  (Fn+l,zn+1)|l  : 

n=m-2  n-m 


C-l 


114  1 

'ttmq  X  62  yields 


V  !||«7."+1||!  +  Atl  ZC  II  •-  K  l  Atll  7,m  II 

1  ~~  b  1 


m-l  ,  2  II  ,  1.2  „  C-l 

♦  l  |||5zn+1|||n  +  l  II  z"*1^  (At)  2  +  |  l  (Fn+l,«zn+1) 

ir-m-2  n-m-2  trm 


(4.7) 


also,  setting  X  =  (n+l)Szn+1  yields 
t-1 


^  (n+l)  HI  6?,n+l  HI  2  ♦  C  Atll  7, Mi;’  <  K(ilmAt  II  7,m  l(  +  Y  II  5  z"  II 


iv  m 


e-i 


n=m 


faAt  ( It  At )  II  n+l 

22  7 


t-1 


n-m-2 


.  +  |  y  (n+l) (Fn+1 ,5zn+1) |l 
n+l  1  ‘ 
a  n  m 


4  .  S) 


Proof :  See  [7]. 

The  following  version  of  the  discrete  Gronwall  lemma  is  ti ivial 


Lemma  4.1.  Let.  f_.  ^  0,  fi .  0,  and  -y  0.  Assume  that  for  n 

n-1 

f  <  T  P.f.At  +  y 
n  —  ,*■  13 


.ii  d 


j*m 

n-1 

l  6. At  <  M  . 

j>*m  1 


Then,  f  y  exp  M,  n  =  m-l,...,t. 


We  shall  assume  that  an  efficient  start-up  procedure  using  the  same  precon¬ 
ditioned  iterative  methods  as  described  in  Section  1  has  been  used  to  determine 
initial  approximations  satisfying 


2 

i-1 


<  K[h 


+  ( At )  M 


(4  .<!) 


Foi  the  description  of  such  a  start-up  procedure  and  proof  of  the  given  estimates, 
,.ee  f  3 1  . 
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We  next  state  the  major  result  of  the  paper. 


Theorem  4.1.  l.et  u  and  U  satisfy  (4.1)  and  (4.2),  respectively.  Let  V  be  t  h, 

iterative  variant  of  U  satisfying  (4.^),  (3.10),  and  (3.11)  with  p  satisfying 

,  r  A  1  n 

(4.21)  below.  Let  u  <  L  (J»Hr)  n  W4(J(W,)  and  either 


*)  /  t|g^(’,t)  ‘dt  <_  K  when  0  is  H*-regular  and  h‘  <_  CAt , 


or 


0 


3  vi  2  i* 

b)  c  L  (JjH  ) 


Then  there  exist  constants  K^(u),  depending  upon  the  norms  of  u,  and  h  (  and 


r0  such  that  if  r  >  d/2,  At  <  min{i0,hd/6} ,  and  h  <  h0. 


»upltun  -  vnll  2i^8(u)lhr  4  (At)  '*1 
n 


.  n+1 


Proof :  Letting  X  -  C  in  (4.3)  with  m  -  3  and  using  (4.n),  we  obtain 

~\lUtn+l  It  4Atllrn+lH' 

n+l 


1^4  l  At»Cn+lfi"  ,) 

n-3  an+1 

<  K4  lilt.3  II2  4  Y<lkn+1i2  4  II  cnM  f)  ♦  lY  l  Tfl(Cn+1)|]  . 


(4.10) 


n-3  i-1 


Next,  we  see  that  from  (2.7)  and  (2.11), 


a «  llc-1  * 


(4.11) 


«-l. 


<  Kq(u)  h2r  4  I  l  IUn+1|fn+1At  , 


n-3 


where 


s  ‘  *  IHwV 


We  note  that  use  of  (2.7b)  instead 


of  Lemma  (2.2),  would  have  required  the  assumption  e  L*(JjHr),  a  much  stronger 

dt 

smoothness  assumption.  From  Laima  (4.1)  we  see  that 


n"J  n-3  a 


(4.1  ) 


We  next  use  (2.4)  and  smoothness  of  W  to  obtain  the  bound 


12- 


] 


e-i 


l  |T"*1R"*l)|  <  K  l  (U"*l|  *  K,(At)5  *  I 


n=  3 


n=3 


j=0 


(4.13) 


<  K(u)  {  (At) b  +  l  [IUn+lH  +  II  6^n+1|l  lAt)  +  l  l!cn+1l!  At  . 
—  ,  1  8  *•„  n+1 

n=l 


8-1  i 

I  I 

n=3 


Using  (3.8),  (3.11)  and  (3.12)  we  see  that 

'iVl«“l>i  <  Tniv"*1  -vn4ii 

n=3 


n=3 


n+1 


64Vn+1| 


e-i 

£  l  p 

n=3 

t-1 

£  l  K(u)p>  {  l  HI 6C 
n=3  i=0 


n+1 1 


n+1 1 


n+l-i I 


(4.14) 


U  *  <«>4>lll<"*lllln 


£  K  (u)  [  (At)  6  +  l  II  Cn  II  At]  +  i  l  II  Cn  I  n+1At 
n=3  n=3  a 


l-l  p1  , 
♦  l  n+1 


n=0 


16At 


fit 


n+1 1 


Noting  that  the  multiplier  in  the  last  term  on  the  right  side  of  (4.14)  is  bounded 
by  (n+l)/16  usinq  (3.12),  we  combine  (4.10)  -  (4.14)  and  use  (4.9)  to  obtain 

ILM|2  +  i  Y[lkn+1f  +  l5n+1«2  At] 

J  n+1 


n=3 


<  K (u)  [h^r  +  (At)  b  +  l  lkn  II  At] 

n=3 

*  %  Yiii«cn,1iiL  *  T  ^  ii 

n=3  n=3 


(4.15) 


6C 


n+1 1 


We  note  that  if  we  can  bound  the  last  two  terms  on  the  right  of  (4.15),  we  can  then 
use  the  discrete  Gronwall  Lemma  to  obtain  our  result.  In  order  to  bound  the  next 

to  the  last  term  on  the  right  side  of  (4.15)  we  let  x  =  6j;n+1  in  (4.3)  and  use 
(4.7)  to  obtain 


8-1 

I 

n=3 


I  lll^n+1||ln  +  <  K5  lAtll  c3  llj  +  |!|6C2|||1  +  III  6c' 


(4.16) 


8.-1  4 


+  YllWcAt)2  +  |Y  l  T"+1(6Cn+1)|] 
n=l  n=3  i*l 
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As  in  ( 4 . 1 1 >  we  use  (2.7)  and  (2.11)  to  obtain 


l-l 


,  ,4-1  ,2 

2c  l  r  in  ..n+l 


l  l-lfWM  <  K(u)h^r  +  i  l  III  6c 
n« 3  n- 3 


(4.17) 


Similarly  we  see  that 

<  KluHltl6  *i*il|ll«tWl|ll’  ' 
n*3  n=3 

Using  (2.4)  we  then  see  that 

tiVfl(«CB+l>|  <  *  T<Un+1|  +  K(At)3  +  l  |«cn+1*j|)|«cn+1|At 


(4. IS) 


n«3 


n=3 


j*0 


(4.19) 


^'i1  Ilk"*1  III!  ♦  Kiuuut)6 .  'i‘iir'r« . 

n=l  n=l 

Then,  as  (4.14),  we  use  (3.8),  (3.11),  and  (3.12)  to  obtain 

n+1  m 


T |T;+1(6cn+1)|  <  ll  P;|||iV+1| 

n=3  n=3 


lnlH4«  "■  n 


4-1 


1  l  “lO.n^ 


n=3 


i=0 


+  1'iHln_i  +  (At)4}  HI  6^n+1  |||n  (4.20) 


4-1 


<  K (u) (At)  +  12  y  p'K 

—  t-  r  n 


n=0 


n  10,n 


n+1 


where  K. „  depends  upon  local  upper  and  lower  bounds  for  the  coefficients  a 
10, n  n 

and  c  (see  (4.21)).  Then  iterating  on  the  preconditioned  iterative  procedure 
n 

sufficiently  often  that 


0n  <  <49Kin.n> 


-1  =  min{a(t3)  ;  j=n+l ,n,n-l ,n-2 ) 


30,n  48  sup{a(t^)  :  j=n+l ,n , n-1 ,n-2 } 


(4.21) 


combining  (4.16)  -  (4.20),  and  using  (4.9)  we  see  that 
4-1  in  . .  m  2  „  .  „2 


n=3 


Okn*1llln  *  Ati.'ii, 

<  K(u)lh2r  +  (At)6  +  l  {(Ln+1|  At  +  IUn+1||  (At)2}] 

n=3 


(4.22) 


In  order  to  bound  the  last  term  on  the  right  side  of  (4.15),  we  let  X  =  (n+l)6i; 
and  use  (4.8)  to  obtain 


n+1 
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T<„.n||k"*‘ir  .  ■  K.MAtIU'f  .  'j‘lk" 

n-1  n  1  1  n-1 


f-1 

*  l 


1' 
n-  1  a 


f-1  4 


♦  I  l  l  T^Ntn.Dfirr1) 


(4.21) 
'>♦1, 


n  3  i<=l 


Wf  not  i  that  (4.  ')  and  (4.2  )  can  be  used  to  Ixuind  t  he  first  term  on  the  r  ight  side 
(  !  (4.2  i)  .  We  next  obtain 


f-l 

i  V  ,  . , _n+ 1  p+1 

i  l  IntilT,  (  ‘ 

n-3 


1 


e-i 

l 

n  \ 


)  I  i  K  l  lllnn+lll  IUcnflll  ♦ 


l  •|dtn 


i=o 


n+1‘jll_1IUf.n+1H1)  (ntl)At 


(4.24) 


1  l  (n+1 )  III  Af,n+l  |l  +  k  Y  II  nn+ 1  II  At  +  K  £  Id  nn+ 1  II  (n+l)At  . 

n-3  n  n-3  n=l  * 


If  n  ■  L>  (J;H  ),  we  have  from  (2.7)  that 


Y«nn+l“2 

n=  1 


At  <  K (u) h 


2r 


(4.2S) 


Then,  usinq  (2.11)  we  have,  if  h~  <  At , 


t  *  TrlW .  b"*‘|L 

n-1  n=l 


2r+2 


(n+1) At 


<  k </  t in u ( - , t )  ii2  +  ||^c,t)  ||;:idt)h2r 


3t 


(4.26) 


Note  that  h  <_  At  is  not  a  stronq  restriction  for  these  hiqh  order  time-steppinq 
methods .  The  constant  on  the  right  of  (4.26)  determines  the  smoothness  assumptions 
,  3u 

wo  need  on  u  and  7^7  for  this  arqument .  We  note  that  for  linear,  t ime -dependent 
problems  the  assumption 


/  t|||^-(  •  ,t)  IPdt  <  K 
0  r 


(4.27) 


» roughly  equivalent  to  —  c  L  (J;Hr  S,  the  assumption  needed  for  (4.11),  and 


wi.  h  weaker  than  the  assumption  ^  c  L2(J;ll')  which  has  been  made  in  [7,  11,  17] 
f  ,r  similar  estimates.  Usinq  (4.4)  we  see  that 


f-1 


e-i 


]>  |T"+1((n+l)5(;n+1)  I  <  ^  Y  (n+1)  6tn+1  „  +  K(u) (At) 6  .  (4.28) 


n=3 


n=3 
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We  next  consider  the  T3  term  from  (4.23).  Note  that 

iY  Tf‘un.usi"*iii  <  iY<f  sv*i,«c"*i)(„.1)iti 

n=3  n=3 

*  lY<^  «V*1.«c"*1>(n.l>M:l  *  |Y<|a  rl.«r1xn.l,4t! 

n»  3  n*  3  u 

-  T5  +  T6  +  T7  * 

can  be  bounded,  using  (2.4),  as  follows 

T,  1  K  l  l  |6cn+1"j|  |«Cn+1| (n+l)At 
n=3  j=0 

1  T7-  l  l«Cn+1H  (n+1)  +  K  X  I  «Cn+1  H  (n+1)  (At) 2  . 

_ 1  _ «  A 


(4.29) 


(4.30) 


Then  using  a  technical  summation  by  parts  argument  and  estimates  like  those  used  in 
(4.  30  we  can  obtain  (see  112,  p.  27-29)  for  details) 

t5  +  t7  -Te  1  {^n+1H  (n+1>  +  «cn+1«  n+1At>  +  jjIUMI^At 

n»l  a 

+  Kile3  I  At  +  (At) 6  +  l  ill  «cn  1  jn  +  At]  llcn+1ll  +  ft  c”*1  f  (n+1)  At]*  ’  ^ ' 

nal  T 


+  K12ll;Aft  lAt  . 


As  in  (4.20),  we  see  that 

l  |TJ+1((ntl)6cn+1)|  <  I  K10 l  |||«cn+1’i||ln.i  +  (At)4}|||6;n+1|||n(n+l) 

n-3  n=3  i*0 


i  K(u)  <4t)6  *  4Kl0  n  Y^o.Jll^llY-ll 

n*o 


(4.32; 


Next,  by  iterating  sufficiently  often  to  satisfy  (4.21),  combining  (4.23)  -  (4.32), 
and  using  (4.9),  we  obtain 
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n=  1  ’ 


*-lM 

III 

n=3 


+  K(u)  [h2r  +  (At)6]  +  K12lcMI  fAt 


(4.33) 


+  K  (u)  l  {ILn+MI  +  li;n+1|l  (n+1)  At)(H  5;nll  an  ♦  At)  . 

n=l  1 


Now  adding  inequalities  (4.15)  and  (4.33)  to  times  inequality  (4.22)  and  sim¬ 

plifying  we  obtain 


cT  *  IU*l(*At  +  Tnil«Cn+1HNn+l)  +  .  At } 

1  n  n+1 

n=3 


<  i<(u){h2r  +  (At)6}  +  4K  IUMI  *At 


(4.34) 


+  K(u)  Y{IUn+1»2  +  Hcn+1fl'nAt}{«6cn||2n  +  At}  . 

“  I  oo 

n=l 


We  next  indicate  how  to  treat  the  term  multiplied  by  4K  _  on  the  right  side  of 
(4.34).  Note  that  for  some  >  0, 


I  C°+1  fi2  (n+1)  At  -  IUn«2nAt  -  ft  ^  l\t 

<  e1(n+l)ll5cn+1ft  +  kJI  c"  I  At  . 


(4.35) 


We  sum  (4.35)  from  n=3  to  n  =  l  -  1 ,  multiply  the  results  by  4K  +  and 

-1  A  2  2 

add  the  final  inequality  to  (4.34).  Then  take  e  <_  (8K  +  1)  .  Next,  we  make 

the  induction  hypothesis  that 


t-i,.  .,2 

l  II  6;n  II  Bn  <  1  . 

n=l  L~ 


(4.36) 


Then  it  follows  from  (4.34)  -  (4.36)  and  Lemma  4.3  that 

£-1..  ..2 

l  II  6cn  II  n  £  2  exp{  (l+T)K(u)  }K(u)  [h2r  +  (At)6] 
n=l 


(4.37) 


It  then  follows  from  (4.37)  and  the  inverse  hypothesis  (2.3.c)  that 

I  H  <5^"  I  <  K  h'd  l  II  6^n  II  n  £  Kh'd[h2r  +  (At)6]  .  (4.38) 

n=l  L  n=l 

We  note  that  the  right  hand  side  of  (4.38)  tends  to  zero  as  h  tends  to  zero  if 
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r 


d 


i-j  .  ;*•' 


—  and  At 


d 

t> 


which  lust  it  ms  the  induction  hypothesis.  Since  this  implies 

He1  i  <_  K|h*'  ♦  (At)*']  ,  (4.40) 

the  result  follows  from  (4.1  ),  Uwru  2.1,  .tiui  the  triangle  i  nequa  1  i  t  y . 

We  note  that  similai  theorems  hold  for  the  original  nonlinear  problem  and  for 
the  other  various  multistep  methods  presented.  Also,  if  i]  is  a  rectangle,  rec¬ 
tangular  solid,  or  unions  of  these  regions,  alternating  direction  variants  ot  the 
multistep  methods  presented  here  are  even  more  computationally  efficient.  See 
(4,  5]  for  these  results. 

V.  ('omput  at  i  on.i  1  dons  idi;i  at  ion-- .  In  this  sect  ion  we  shall  consider  some  rough 
operation  counts  to  estimate  the  computat ional  complexity  of  the  methods  presented 
here.  We  shall  see  that  the  preconditioned  iterative  methods  allow  us  to  obtain 
very  nearly  optimal  order  work  estimates  and  are  thus  very  efficient  computation¬ 
ally  . 


We  shall  give  estimates  for  d  *  2.  The  procedures  ot  setting  up  and  factoring 
n  .  3/  2 

I,  requires  0(M  )  operations,  where  ^  dim  Jt.  The  solution  of  (3.2),  given 

the  factorization,  requires  0(M  log  M)  operations.  Such  bounds  have  been  shown 
to  be  minimal.  If  we  conjecture  the  validity  of  t he  above  estimates  for  our  problem 

and  refactor  L  and  solve  (3.2)  at  each  time  step,  the  total  amount  of  work,  done 
is 

3/2  3/2 

0(N{M  -*■  M  loq  M})  -  0(NM  )  ,  (5.1) 

where  N  is  the  total  number  of  time  steps  (N  *  (At)  S.  Note  that  the  work  of 
factorization  dominates  the  estimate. 

Using  the  preconditioned  iterative  procedure  presented  here,  only  the  precon¬ 
ditioner,  L  ,  must  be  factored.  Let  ic  be  the  number  of  iterations  needed  to 

n 

achieve  the  necessary  norm  reductions  in  ( 3 . 1 1  >  and  (3.12).  We  note  that  k  can 

n 

be  bounded  by  a  fixed  constant  e  which  is  independent  of  h,  n,  and  At.  Using 
this  method  the  total  work  done  is 

0(m’  2  ♦  NkM  log  M)  .  (5.2) 

Since  balancing  the  spatial  and  temporal  errors  yields 

r 

N  *  (At)*1  *  h  11  •  0(M*  “') 

we  note  that  for  r  ^  y ,  the  work  of  solving  dominates  the  estimate,  while  for 
r  <  u  the  amount  of  work  of  solving  is  even  less  than  the  work  to  factor  one 
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matiix,  a  necessary  piece  of  work,  clearly,  in  any  case,  (5.2)  is  much  pi.-t.  table 
to  (S.l)  .  Also,  since  the  total  numbet  of  unknowns  in  the  pioblem  is 

O ( NM l  , 

is.2)  represents  a  nearly  optimal  ordet  w\>:  k  est  uiutf  when  the  woi  k  is  at  least  a- 
muc*  is  factoring  one  matiix.  It  alternating  diieotion  vat  iant  s  ot  these  met  h.vh 
ca  be  used,  the  log  M  term  can  he  t  emoved  from  (5.2)  and  opt  imal  ordei  work  t  i 
mates  are  obtained  (see  [4,  *M  ). 


It  1 s  computat lonal ly  wasteful  to  iterate  exactly  c  t imes  at  each  t ime  step 

m  oidet  to  achieve  the  pessimistic  bounds  on  p  given  m  (4. .Ml.  Instead,  one 

n 

.in  monitor  the  norm  reduction  actually  produced  at  each  t  ime  step  ot  the  iteiatioi 
>ud  top  iterating  when  sufficient  norm  reduction  is  achieved.  Additional  stopping 
criteria  can  he  imposed  in  this  monitoring  process.  See  I'M  for  a  discussion  ot 
stopping  criteria  for  related  methods. 
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